Sp = [0,1;0,0];
Sm = [0,0;1,0];
I = eye(2);
Z = [1,0;0,-1];

disp(kron(Sp,Sm) + kron(Sm,Sp)+1/2*(kron(Z,Z)+kron(I,I)));
disp(Sp*Sm - Sm*Sp);

S1p = kron(Sp,Z);
S1m = kron(Sm,Z);
S2p = kron(I,Sp);
S2m = kron(I,Sm);
S1z = kron(Z,I);
S2z = kron(I,Z);

disp(-S1p*S2m + S1m*S2p + 1/2*(S1z*S2z + kron(I,I)));

c1 = S1m;
c2 = S2m;

disp(-c1'*c2+c1*c2'+1/2*((c1'*c1 - c1*c1')*(c2'*c2 - c2*c2') + kron(I,I)));
disp(-c1'*c2-c2'*c1-c1'*c1-c2'*c2+eye(4)+2*c1'*c1*c2'*c2);


T = [1,0,0,0;0,0,1,0;0,1,0,0;0,0,0,-1];

disp(T'*c1*T - c2);
disp(T'*c2'*T - c1');

disp(-c1'*c1-c2'*c2-c1'*c2-c2'*c1 + eye(4));


Num = 8;
c = fn_cc(Num);
cp = fn_ccp(Num);

Pf = fn_Pf(c);
vacuum = fn_vacuum(Num);

Kp = 1;
for i = 1:(Num-1)
    Kp = Kp*Pf{i,i+1};
end
Kp = (c{1} + c{1}')*Kp;
K = Kp';

for i = 1:(Num-1)
    fprintf('i = %d, error = %e \n', i, norm(Kp*c{i}*K - c{i+1}));
end
fprintf('i = %d, error = %e \n', Num, norm(Kp*c{Num}*K + c{1}'));

fprintf('the invariance of the vacuum %e \n',...
    norm(Kp*vacuum - (c{1} + c{1}')*vacuum));

bqd = cell(1,Num);
bq = cell(1,Num);
bqdp = cell(1,Num);
bqp = cell(1,Num);

for q = 1:Num
    
   tmp = 0;
   
   for i = 1:(Num-q)
       tmp = tmp + c{i+q}'*c{i};
   end
   
   bqd{q} = tmp;
   
   for i = 1:q
       tmp = tmp + c{i}'*c{Num-q+i};
   end
   
   bqdp{q} = tmp;
   
   tmp = 0;
   
   for i = 1:(Num-q)
       tmp = tmp + c{i}'*c{i+q};
   end
   
   bq{q} = tmp;
   
   for i = 1:q
       tmp = tmp + c{Num-q+i}'*c{i};
   end
   
   bqp{q} = tmp;
    
end

fprintf('test bq : \n');
for q = 1:(Num-1)
    fprintf('q = %d, error = %e \n', q, ...
        norm(Kp*bq{q} - bq{q}*Kp + (c{1}'*c{1+q} ...
        + c{Num-q+1}'*c{1}')*Kp));
end


disp(norm(bq{1}'*bq{2}-bq{2}*bq{1}'));
disp(norm(bq{1}*bq{3} - bq{3}*bq{1}));
disp(norm(bq{2}*bq{3} - bq{3}*bq{2}));
disp(norm(bq{2}'*bq{3} - bq{3}'*bq{2}));

fprintf('the important: the rank of the boson operator. \n');

for q = 1:Num-1
    
   m = bq{q};
   
   [v,d] = eig(m);
   
   fprintf('Num = %d, q = %d, rank of vector = %d \n', ...
       Num, q, rank(v));
    
end

tmp = bq{2}*bq{1}' - bq{1}'*bq{2} ...
    - c{1}'*c{2} + c{7}'*c{8} ;

disp(norm(tmp));

q1 = 2;
q2 = 3;

LHS = bq{q1}*bq{q2}' - bq{q2}'*bq{q1};
RHS = 0;

for i = 1:q1
    RHS = RHS + c{q2-q1+i}'*c{i};
end

for i = 1:q1
    RHS = RHS - c{Num-q1+i}'*c{Num-q2+i};
end

disp(norm(RHS - LHS));


psi = cell(1,Num);

L = 1.2;

for ix = 1:Num
    
   tmp = 0;
   for ik = 1:Num
       tmp = tmp + 2*pi/L*exp(-1j*2*pi/L*ik*L/Num*ix)*c{ik};
   end
   
   psi{ix} = tmp;
    
end


for q = 1:Num
    for x = 1:Num
        tmp = norm(bqp{q}*psi{x}-psi{x}*bqp{q} ...
            +exp(1j*2*pi/L*q*L/Num*x)*psi{x});
        fprintf('x = %d, q = %d, error = %e \n', x, q, tmp);
        
        LHS = bq{q}*psi{x}-psi{x}*bq{q} ...
            +exp(1j*2*pi/L*q*L/Num*x)*psi{x};
        RHS = 0;
        for n = 1:q
            RHS = RHS + 2*pi/L*c{n}*exp(-1j*2*pi/L*(n-q)*L/Num*x);
        end
        
        fprintf('check the formula with the exact bq : %e \n', ...
            norm(LHS - RHS));
        
    end
end


vx0 = [1,1,1,1,0,0,0,0];
N = length(vx0);
nx = 1;
L = 1;
psi = fn_spin_wave_base(nx,L,vx0);
disp(psi.amp);
disp(psi.base);

vx0 = [1,1,0,1,0,0,1,0,1,1];
psi = fn_left(vx0,2);
disp(psi.amp);
disp(psi.base);

L = 1.2;
vx0 = [1,1,1,1,0,0,0,0];
N = length(vx0);
fprintf('check the exact formula: \n');
for nx = 1:N
    for nq = 1:N
        
        psi_0.amp = 1;
        psi_0.base = vx0;
        
        psi_1 = fn_spin_wave(nx, L, psi_0);
        psi_1 = fn_annihilation(psi_1,nq);
        psi_2 = fn_annihilation(psi_0,nq);
        psi_2 = fn_spin_wave(nx, L, psi_2);
        
        psi_2.amp = - psi_2.amp;
        LHS = fn_add(psi_1,psi_2);
        
        RHS1 = fn_spin_wave(nx, L, psi_0);
        RHS1.amp = -exp(1j*2*pi/N*nq*nx)*RHS1.amp;
        RHS2.amp = [];
        RHS2.base = [];
        
        for n = 1:nq
            
            tmp = fn_ani(psi_0, n);
            tmp.amp = 2*pi/L*exp(-1j*2*pi/N*(n-nq)*nx)*tmp.amp;
            RHS2 = fn_add(RHS2,tmp);
        end
        
        RHS = fn_add(RHS1, RHS2);
        
        
        fprintf(' nq = %d, nx = %d  \n', nq, nx);
        fprintf(' error = %e \n', norm(RHS.base - LHS.base) ...
            + norm(RHS.amp - LHS.amp));
        fprintf('============================\n');
    end
end

psi_0.amp = 1;
psi_0.base = vx0;

disp(fn_qn1(psi_0));
disp(fn_qn2(psi_0));

psi_1 = fn_creation(psi_0,2);
disp(fn_qn1(psi_1));
disp(fn_qn2(psi_1));

psi_2 = fn_annihilation(psi_1,2);
disp(fn_qn1(psi_2));
disp(fn_qn2(psi_2));

vx = floor(2*rand(1,Num));
psi.amp = 1;
psi.base = vx;

disp(psi.base);

for kk = 1:40
    fprintf('kk = %d \n', kk);
    psi = fn_annihilation(psi,2);
    disp(fn_qn2(psi));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function psi = fn_creation_ntimes(psi_0,nq,n)

psi = psi_0;

for i = 1:n
    psi = fn_creation(psi,nq);
end

end

function res = fn_qn1(psi)
    vx = psi.base;
    res = sum(vx,2);
end

function res = fn_qn2(psi)
vx = psi.base;
L = size(vx,2);
M = size(vx,1);

tm = ones(M,1)*(1:L);
res = sum(tm.*vx,2);
end

function psi = fn_spin_wave(nx, L, psi0)

psi.amp = [];
psi.base = [];

amp = psi0.amp;
base = psi0.base;

for n = 1:length(amp)
    
    tmp = fn_spin_wave_base(nx, L, base(n,:));
    tmp.amp = tmp.amp*amp(n);
    psi = fn_add(psi, tmp);
    
end

end


function psi = fn_spin_wave_base(nx,L,vx0)

psi.amp = [];
psi.base = [];
N = length(vx0);

for nk = 1:N
    
    psi_single = fn_single_ani(vx0,nk);
    psi_single.amp = psi_single.amp*2*pi/L*exp(-1j*2*pi/N*nk*nx);
    psi = fn_add(psi,psi_single);

end

end


function psi = fn_ani(psi0,n)

amp = psi0.amp;
base = psi0.base;

Num = length(amp);

psi.amp = [];
psi.base = [];

for i = 1:Num
    
    tmp = fn_single_ani(base(i,:),n);
    tmp.amp = tmp.amp*amp(i);
    
    psi = fn_add(psi, tmp);
    
end

end


function psi = fn_single_ani(vx,n)

if vx(n) == 0
    psi.amp = 0;
    psi.base = vx;
    return;
end

vx(n) = 0;
amp = (-1)^sum(vx((n+1):end));

psi.amp = amp;
psi.base = vx;

end


function res = fn_Pf(c)

Num = length(c);

N = size(c{1},1);
I = eye(N);
res = cell(Num,Num);

for i = 1:Num
    for j = 1:Num
        if i==j
            res{i,j} = I;
        else
            c1 = c{i};
            c2 = c{j};
            res{i,j} = I - c1'*c1 - c2'*c2 - c1'*c2 - c2'*c1;
        end
    end
end

end

function res = fn_YN(Num)

tv = (-1).^(0:(Num-2));
m = diag(tv,1) - diag(tv,-1);
m(Num,1) = m(Num,1)-1;
m(1,Num) = m(1,Num) + 1;

res = -1j*m;

end

function res = fn_ZN(Num)

res = diag((-1).^(0:(Num-1)));

end

function res = fn_XN(Num)

res = [zeros(Num/2), eye(Num/2); eye(Num/2), zeros(Num/2)];

end




function res = fn_C(eta, Neta, cc, ccp)

Num = size(cc,1);
M = size(cc,2);

if Neta == 0
    res = eye(2^(Num*M));
elseif Neta > 0
    res = 1;
    for i = 1:Neta
        res = ccp{i + Num/2, eta}*res;
    end
else
    res = 1;
    for i = -1:-1:Neta
        res = cc{i + Num/2 + 1, eta}*res;
    end
end

end


function [mat, mat_test] = fn_test_commutation_rule(cc, ccp)

mat = cell(length(cc), length(ccp));

for i1 = 1:length(cc)
    for i2 = 1:length(ccp)
        
        c = cc{i1};
        cp = ccp{i2};
        mat{i1, i2} = c*cp + cp*c;
        
    end
end

mat_test = zeros(size(mat));

for i1 = 1:length(cc)
    for i2 = 1:length(ccp)
        
        ma = mat{i1, i2};
        if i1 == i2
            mat_test(i1,i2) = norm(ma - eye(size(ma,1)));
        else
            mat_test(i1,i2) = norm(ma);
        end
        
    end
end
        
end


function res = fn_vacuum(Num)

v = [0;1];

res = 1;
for i = 1:Num
    res = kron(res, v);
end

end


function cc = fn_cc(Num)

cc = cell(1, Num);

for i = 1:Num
    cc{i} = fn_M(i, Num);
end

end

function ccp = fn_ccp(Num)

ccp = cell(1, Num);

for i = 1:Num
    ccp{i} = fn_P(i, Num);
end

end

function res = fn_P(ind, Num)

P = [0,1;0,0];

res = fn_op_jw(P, ind, Num);

end


function res = fn_M(ind, Num)

M = [0,0;1,0];

res = fn_op_jw(M, ind, Num);

end


function res = fn_op(op, ind, Num)

mI = eye(2);

res = 1;

for i = 1:Num
    
    if i == ind
        res = kron(res, op);
    else
        res = kron(res, mI);
    end
    
end

end


function res = fn_op_jw(op, ind, Num)

mI = eye(2);
Z = [1,0;0,-1];

res = 1;
for i = 1:Num
    
   if i<ind
       res = kron(res,mI);
   elseif i == ind
       res = kron(res, op);
   else
       res = kron(res, Z);
   end
    
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function res = fn_base2ind(vx)

tv = 2.^(0:(length(vx)-1));
res = sum(vx.*tv);

end


function psi = fn_apply_to_vacuum(vn, Num)

vx0 = [ones(1,Num), zeros(1,Num)];
psi0.amp = 1;
psi0.base = vx0;

psi = psi0;

for i = 1:length(vn)
    
   n = vn(i);
   if n > 0
       psi = fn_bdmn(psi, i, n);
   end
    
end

end

function res = fn_find_number(M, N)

if M == 1
    res = N;
    return;
end

res = zeros(2^M, M);
pt = 0;
for i = 0:floor(N/M)
    
   res_red = fn_find_number(M - 1, N - M*i);
   if ~isempty(res_red)
       tmp = [res_red, i*ones(size(res_red,1),1)];
       pt = pt + size(res_red,1);
       res(pt - size(res_red,1) + 1 : pt, :) = tmp;
   end
   
end

res = res(1:pt,:);

end


function res = fn_find_all_base(Num, M)

m0 = (1 + Num)*Num/2;

m = M + m0;

t = nchoosek(1:(2*Num), Num);

tsum = sum(t, 2);

tm = t(tsum == m, :);

res = zeros(size(tm,1), 2*Num);

for i = 1:size(tm,1)
   
    res(i,tm(i,:)) = 1;
    
end

end


function psi = fn_bdmn(psi0, m, n)

psi = psi0;
for i = 1:n
    psi = fn_creation(psi, m); 
end

end


function psi = fn_creation(psi0, n)

psi.amp = [];
psi.base = [];

amp = psi0.amp;
base = psi0.base;

L = length(amp);

for i = 1:L
    
   vx = base(i,:);
   phi = fn_right(vx, n);
   phi.amp = amp(i)*phi.amp;
   
   psi = fn_add(psi, phi);
    
end

end


function psi = fn_annihilation(psi0,n)

psi.amp = [];
psi.base = [];
amp = psi0.amp;
base = psi0.base;

L = length(amp);

for i = 1:L
    
    vx = base(i,:);
    phi = fn_left(vx, n);
    phi.amp = amp(i)*phi.amp;
    
    psi = fn_add(psi, phi);
end

end


function phi = fn_right(vx, n)

L = length(vx);
res = zeros(L, L);
amp = zeros(1, L);

pt = 0;

for i = 1:L
    if i + n <=L
       if  vx(i) == 1 && vx(i + n) == 0
          tmp = vx;
          tmp(i) = 0;
          tmp(i + n) = 1;
          num_of_one = sum(vx((i+1):L)) + sum(vx((i+n+1):L));
          
          pt = pt + 1;
          res(pt,:) = tmp;
          amp(pt) = (-1)^num_of_one;      
       end
    end
end

res = res(1:pt,:);
amp = amp(1:pt);

phi.amp = amp;
phi.base = res;

end


function phi = fn_left(vx, n)

L = length(vx);
res = zeros(L,L);
amp = zeros(1,L);

pt = 0;

for i = 1:L
    if i+n<=L
        if vx(i) == 0 && vx(i + n) == 1
            tmp = vx;
            tmp(i) = 1;
            tmp(i + n) = 0;
            num_of_one = sum(vx((i+1):L)) + sum(vx((i+n+1):L));

            pt = pt + 1;
            res(pt,:) = tmp;
            amp(pt) = (-1)^(num_of_one + 1);
        end
    end
end

res = res(1:pt,:);
amp = amp(1:pt);

phi.amp = amp;
phi.base = res;

end


function psi = fn_add(psi1, psi2)

amp = [psi1.amp, psi2.amp];
base = [psi1.base; psi2.base];

psi.amp = amp;
psi.base = base;
psi = fn_combine(psi);

end


function psi = fn_combine(psi1)

amp = psi1.amp;
base = psi1.base;

Num = length(amp);

for i = 1:(Num - 1)  
    if amp(i)~=0
        vec0 = base(i,:);
        for i2 = (i + 1):Num
            vec1 = base(i2,:);
            if max(abs(vec1 - vec0)) == 0
                amp(i) = amp(i) + amp(i2);
                amp(i2) = 0;
            end
        end
    end
    
end

ind = find(abs(amp) < 1e-12);

amp(ind) = [];
base(ind,:) = [];

psi.amp = amp;
psi.base = base;

end